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Abstract 

An  algorithm  is  presented  whieh  uses  the  ehannel  paeket  method  (CPM)  to 
eompute  S-matrix  elements.  The  standard  approach  to  solving  for  the  S-matrix  elements 
is  to  use  states  which  contain  only  positive  or  negative  momentum.  In  doing  so  however, 
the  standard  approach  fails  to  produce  accurate  S-matrix  elements  at  low  energies  since 
accurate  results  are  only  obtainable  over  the  energy  range  defined  by  the  states. 

Therefore,  in  order  to  obtain  accurate  results  at  low  energy  one  must  formulate  states 
which  contain  both  positive  and  negative  momentum.  In  order  to  incorporate  states  which 
have  positive  and  negative  momentum  a  four-by-four  matrix  containing  the  momentum 
expansion  coefficients  of  the  states  is  introduced.  The  approach  does  not  consider 
scattering  from  one  side  or  the  other,  rather  it  considers  both  incoming  and  outgoing 
states  from  the  left  and  right  simultaneously.  Therefore,  during  one  simulation  all  four  S- 
matrix  elements,  and  S_^  are  computed.  Numerical  simulations  of 

the  algorithm  are  carried  out  on  a  conventional  desktop  computer  and  compared  to  the 
analytic  solution  of  the  transmission  and  reflection  functions  for  a  square  and  trapezoidal 
well.  The  simulated  results  agree  very  well  with  the  known  solution  up  until  very  low 
energies,  after  which  the  results  begin  to  oscillate  about  the  theoretical  values.  The 
explanation  of  the  oscillations  is  presented  as  well  as  two  error  analyses:  the  first  showing 
that  the  oscillations  approach  zero  as  the  reactant  and  product  wave  packets  are 
propagated  away  from  the  interaction  region,  and  the  second  showing  convergence  to  the 
correct  S-matrix  elements  for  a  trapezoidal  well.  The  results  indicate  that  if  the  exact 
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M0ller  states  are  used  the  algorithm  will  produee  the  eorreet  S-matrix  elements  aeross  the 
entire  energy  range. 


V 


Acknowledgments 


I  would  like  to  express  a  sincere  gratitude  to  my  faculty  advisor,  Dr.  David  Weeks, 
for  his  guidance  and  expertise  throughout  the  work  of  this  thesis.  His  enthusiastic 
demeanor  during  our  technical  discussions  made  learning  about  scattering  theory  a  joyful 
experience. 

A  special  thank  you  goes  to  my  girlfriend  for  her  constant  love,  compassion,  and  patience 
throughout  this  thesis.  I  would  like  to  thank  her  for  understanding  that  even  during  visits 
it  may  have  seemed  as  though  we  were  still  apart  during  the  intense  research  phase.  I 
would  also  like  to  thank  my  parents  for  always  believing  in  my  dreams  and  their 
continued  support  for  the  choices  I  make. 


Clint  M.  Zeringue 


Table  of  Contents 


Page 

Abstract . iv 

Acknowledgments . vi 

List  of  Figures . ix 

List  of  Tables . xii 

I  Introduction . 1 

1.1  Overview . 1 

1.2  Organization . 4 

II  Time -Dependent  Scattering  Theory . 6 

2.1  The  S-Matrix . 6 

2.2  The  Channel  Packet  Method . 9 

2.2.1  Moller  Operators . 9 

2.2.2  The  Split  Operator  Method . 15 

2.2.3  Computing  the  S-Matrix  Elements . 18 

III  Application  of  the  Channel  Packet  Method . 20 

3.1  Previous  Technique . 20 

3.2  New  Technique . 23 

3.2.1  Summary  of  Steps  for  New  Approach . 27 

3.2.2  An  Interesting  Implication  of  the  New  Approach . 28 

IV  Choosing  the  Parameters . 29 

V  Properties  of  the  New  Algorithm . 33 

vii 


5.1  Verifying  the  New  Approaeh . 33 

5.2  Low-Energy  Seattering  (Square  Well) . 35 

5.2. 1  Diffieulty  in  Construeting  the  Moller  States . 35 

5.2.2  Approximate  Square  Well  Results . 37 

5.2.3  Convergenee  to  Square  Well  based  on  Ar . 40 

VI  Square  Well  Error  Analysis . 43 

6.1  Methodology . 43 

6.2  Results . 44 

VII  Eow-Energy  Seattering  (Trapezoidal  Well) . 47 

7.1  Solving  Sehrodinger's  Equation . 47 

VIII  Trapezoidal  Well  Error  Analysis . 50 

8.1  Results . 52 

IX  Cone  lus  ion . 56 

Appendix  A . 58 

Appendix  B . 60 

Appendix  C . 61 

Bibliography . 62 

Vita . 63 


viii 


List  of  Figures 


Figure  Page 

1 .  Seattering  from  an  arbitrary  loealized  potential  {V(x)-0  exeept  in  Region  II) . 6 


2.  The  initial  produet  and  reaetant  Gaussian  wave  paekets|'p.^)  and  |4'„,)  plotted  in 

the  eoordinate  representation.  The  square  well  potential  depth  is  sealed  by  a 
faetor  of  10.  The  probability  density  is  shown . 11 

3.  The  initial  produet  and  reaetant  Gaussian  wave  paekets|'p.^)  and  plotted 

in  the  momentum  representation  where  only  +k  is  eonsidered.  The  probability 
density  is  shown . 12 

4.  The  probability  density  of  the  intermediate  Moller  states,  14-/)  and  )  (after 

t  =  30000  a.M.)  plotted  in  the  eoordinate  representation.  The  square  well  potential 
is  also  plotted . 14 

5.  The  Moller  states,  |'r+)  and  I'?-)  plotted  in  the  eoordinate  representation.  The 
square  well  potential  has  been  sealed  by  a  faetor  of  10. The  probability  density  is 


shown . 16 

6.  Initial  produet  and  reaetant  Gaussian  wave  paekets  plotted  in  the  momentum 

representation  showing  overlap.  The  probability  density  is  shown . 22 


7.  Initial  produet  and  reaetant  Gaussian  wave  paekets  \  and  \  plotted  in  the 

”'i/  '"■7/ 

eoordinate  representation.  The  square  well  potential  has  been  sealed  by  a  factor 
of  10.  The  probability  density  is  shown . 26 

8.  Initial  product  and  reactant  Gaussian  wave  packets  \  and  .j,  \  plotted  in  the 

^  out  \  j  ^irt]  I 

2/  ij 

momentum  representation.  The  probability  density  is  shown . 26 

9.  (a)  Approximate  square  well  potential  for  large  Ax  .  (b)  Approximate  potential 

showing  the  convergence  to  an  accurate  square  well  as  Ax  is  decreased . 30 

10.  Real  part  of  the  correlation  functions  for  simulation  1.  The  dashed  function 
represents  Ci(0=C3(0  ,  while  the  solid  line  represents 

C2(0  =  . 33 


IX 


Figure 


Page 


1 1 .  Transmission  (upper  )  and  refleetion  (lower)  eurves  for  the  square  well  of 
depth  F  =  -0.01  a.u  and  length  2a  =  4 .  The  dashed  lines  represent  the 
analytie  solutions  while  the  solid  lines  indieate  the  numerieal 
results . 34 


12.  (a)  Approximate  Moller  states  after  t  =  ±800000  {a.u.)  propagation  time. 
Notiee  that  the  wave  funetions  have  spread  dramatieally  .(b)  Zoomed  in 
portion  of  interaetion  region -2  <x<2 . 


13.  Approximate  square  well  (dashed)  and  square  well  (solid)  for 


V{x)  = 


J-0.01  -2<x<2 

[O  elsewhere 


36 


37 


14.  Transmission  [  ]  and  refleotion[  ]  eurves  for  the  square  well 

of  depth  =  -0.01  au.  and  length  2a  =  4 .  The  analytie  solutions  are  (solid) 
while  the  simulated  results  are  (dashed),  (b)  Zoomed  in  portion  of  low  energy 
region  showing  oseillatory  behavior . 39 

15.  Relative  error  vs.  Energy  for  sample  test  for  square  well  of  depth 

Vg  =  -0.01  a.M.and  width  2a  =  4  .  The  energy  value  of  the  measured  data 

point  is  labeled . 40 

16.  Relative  error  vs.  Ax  for  square  well  of  depth  V (x)  =  -0.01  (a.u)  and  width 
2a  =  4  .  An  approximate  fit  reveals  that  the  error  exponentially 

inereases  as  the  eoordinate  grid  spaeing  is  inereased . 41 

17.  5'^j.  and  eurves  showing  the  oseillations  used  to  perform  the  error 

analysis . 45 

18.  Average  amplitude  of  the  live  oseillations  shown  in  figure  17  vs.  the  time 

for  whieh  the  produet  and  reaetants  were  propagated  away  from  the  interaetion 


region . 46 

19.  Trapezoidal  potential  showing  effeet  of  having  a  large  Ax  .  The 

five  regions  are  labeled  under  their  respeetive  parts . 47 


20.  (a)Plot  of  the  Transmission  eoeffieient  vs.  energy  (a.u.)  for  a  square  well 
(solid  line)  and  trapezoidal  well  (dashed)  with  Ax  =  0.1 .  (b)  Plot  showing 
the  eomparison  of  the  trapezoidal  potential  and  a  square  potential. 


X 


Figure 


Page 


(e)  Relative  error  of  the  trapezoidal  transmission  eoeffieient  and  square  well  vs. 
energy  (a.u.) . 48 

21.  (a)Plot  of  the  Transmission  eoeffieient  vs.  energy  (a.u.)  for  a  square 

well  (solid  line)  and  trapezoidal  well  (dashed)  with  Ax  =  a  (Maehine  Epsilon). 

(b)  Plot  showing  the  eomparison  of  the  trapezoidal  potential  and  a  square 
potential,  (e)  Relative  error  of  the  trapezoidal  transmission  eoeffieient 
and  square  well  vs.  energy  (a.u) . 49 

22.  Plot  of  the  T+R  vs.  Energy  (a.u)  showing  that  it  indeed  equals  one  for 

all  energies . 49 

23.  Trapezoidal  potential  well  showing  diserete  points  along  eurve. 

The  well  is  0.01  atomie  units  deep . 50 

24.  (a)  Transmission  [  ,  s_^  ]  and  refieotion[  ,  s_^  ]  eurves  for  the  trapezoidal 

well.  The  analytie  solutions  are  (solid)  while  the  simulated  results  are 
(dashed),  (b)  Zoomed  in  portion  of  low  energy  region  showing  oseillatory 
behavior . 51 

25.  (solid)  and  Sj.  (dashed)  vs.  Time  along  with  power  fits  showing  oonvergenoe....53 

26.  Sehematie  showing  integrated  residue  (shaded)  inside  interaetion  potential 

(Trapezoidal  Well) . 54 

27.  Eigure  28.  Sj^  and  Sj  along  with  the  integrated  data  showing  the  eorrelation 

between  residue  left  inside  interaetion  region  and  error  assoeiated  with 

ringing . 55 


XI 


List  of  Tables 


Table  Page 

1 .  S-matrix  seleetion  table . 21 

2.  Times  of  propagation  for  Error  Analysis  Two . 53 

3.  Sample  parameters  of  simulation  for  =  ±4 . 61 


xii 


A  NEW  APPLICATION  OF  THE  CHANNEL  PACKET  METHOD  FOR  LOW 
ENERGY  1-D  ELASTIC  SCATTERING 


I  Introduction 

1.1  Overview 

On  the  highway,  we  try  to  avoid  eolliding  with  other  ears.  On  the  pool  table,  we 
tailor  eollisions  to  suit  our  purposes.  In  partiele  aeeelerators,  we  interpret 
eollisions  to  reveal  the  innermost  strueture  of  matter.  On  a  eosmie  seale,  we 
wonder  if  a  eollision  between  an  asteroid  and  Earth  led  to  the  extinetion  of  the 
dinosaurs.  Collisions  are  important  features  of  our  physieal  universe. 

Jay  M.  Pasaehoff,  Physics  (1999) 

Classieal  meehanies  deseribes  the  properties  and  effeets  of  various  eollisions 
through  the  conservation  laws  of  energy  and  momentum.  In  the  microscopic  world 
however,  one  must  adopt  the  theory  of  quantum  mechanics  to  characterize  the 
interactions  between  particles.  In  classical  mechanics  one  only  needs  to  specify  the 
velocities  and  masses  of  the  interacting  constituents  to  completely  characterize  a 
collision;  the  quantum  collision  picture  or  “scattering  theory”,  is  a  much  more  subtle 
affair.  In  quantum  mechanics  velocity  and  position  are  viewed  in  terms  of  probabilities. 
The  purpose  of  collision  theory  in  the  quantum  mechanical  regime  is  then  to  calculate  the 
probability  that,  as  the  result  of  a  collision,  a  particle  will  “scatter,”  or  deviate  from  its 
initial  path.  Scattering  theory  was  developed  to  describe  the  statistical  nature  of 
microscopic  collisions. 
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The  Strahlung  matrix  (S-matrix)  is  the  core  of  scattering  theory.  The  S-matrix  is  a 
unitary  matrix  that  relates  the  final  state  in  the  infinite  future  (out-channel)  and  the  initial 
state  in  the  infinite  past  (in-channel).  The  individual  elements  in  the  S-matrix  are  known 
as  scattering  amplitudes;  squaring  the  modulus  of  the  S-matrix  elements  produces  the 
transmission  and  reflection  probability  coefficients  T  and  R  for  an  arbitrary  potential 
V(x) .  The  channel  packet  method  (CPM)  is  a  time  dependent  approach  to  scattering 
theory  that  produces  the  S-matrix  elements  numerically. 

Dr.  David  E.  Weeks  and  Dr.  David  J.  Tannor  developed  the  CPM  in  the  early 
1990’s.  Tannor  and  Weeks  used  Moller  operators  to  propagate  reactant  and  product  wave 
packets  into  interaction  regions  of  step-like  potential  barriers  and  wells  as  well  as  a 
collinear  reaction  of  atomic  hydrogen  and  dihydrogen.  The  S-matrix  elements  were 
calculated  through  a  Fourier  transform  of  a  time-dependent  correlation  function  between 
the  Moller  states.  Using  the  S-matrix  elements  Weeks  and  Tannor  were  able  to  produce 
transmission  and  reflection  curves  that  agreed  well  with  the  known  solutions.  However, 
at  low  energies  the  results  of  their  calculations  deviated  significantly  from  the  theoretical 
predictions.  [1-3] 

A  more  complete  algorithm,  capable  of  yielding  good  results  along  the  entire 
energy  spectrum,  is  needed.  The  goal  of  this  project  is  to  explore  1-D,  elastic  scattering 
and  apply  a  new  approach  to  the  CPM  in  an  attempt  to  correct  the  errors  that  occur  at  low 
energies.  A  new  algorithm  is  developed  to  approach  the  problem  from  a  different 
perspective.  Previous  algorithms  described  in  [1-3]  contain  a  specific  requirement-  the 
reactant  and  product  wave  functions  are  formulated  such  that  they  have  either  all  positive 
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or  negative  momentum.  Weeks  and  Tannor  showed  that  aeeurate  results  were  only 
attainable  over  the  range  of  momentum  for  which  the  amplitude  of  the  product  and 
reactant  wave  packets  is  nonzero.  Since  previous  methods  confined  the  momentum 
amplitude  of  the  wave  packets  to  go  to  zero  near  low  momentum  (and  hence  low  energy) 
information  was  lost  and  as  a  result,  deviations  from  the  analytic  solutions  were  observed 
in  the  low  energy  portions  of  the  transition  curves. 

The  objective  of  this  thesis  is  to  develop  a  new  algorithm  capable  of  producing 
better  results  in  the  low-energy  limit.  To  achieve  this  objective  we  employ  several 
adaptations  to  previous  applications  of  the  CPM.  First,  the  product  and  reactant  wave 
packets  are  constructed  to  contain  both  positive  and  negative  momentum.  As  an  effect  of 
the  previous  requirement,  scattering  from  the  left  and  right  is  considered  simultaneously, 
and  the  entire  S-matrix  is  produced  during  a  single  propagation.  Lastly,  a  significant 
increase  in  simulation  time  as  well  as  grid  size  is  required  to  propagate  the  product  and 
reactant  wave  packets  into  the  asymptotic  channels.  The  last  requirement  stems  from  the 
first  in  that  with  both  positive  and  negative  momentum  propagation  out  of  the  interaction 
region  becomes  more  time  consuming  do  to  spreading  of  the  wave  functions  in  both 
directions.  The  nature  of  the  approach  taken  in  this  thesis  requires  that  the  wave  packets 
contain  both  positive  and  negative  momentum.  To  truly  propagate  the  wave  packets 
completely  out  of  the  interaction  region  requires  an  infinite  grid.  Since  this  requirement 
cannot  be  met  one  must  accept  this  limitation  and  attempt  to  propagate  far  enough  so  that 
the  amplitude  remaining  inside  the  interaction  region  is  negligible. 
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There  are  two  main  contributions  presented  in  this  thesis:  an  algorithm  that 
produces  a  significant  improvement  in  the  reflection  and  transmission  curves  at  low 
energies,  and  a  program  developed  to  analytically  solve  for  the  transmission  and 
reflection  coefficients  for  a  trapezoidal  well.  In  addition,  the  explanations  of  problematic 
results  at  low  energies  are  characterized. 

1 .2  Organization 

This  thesis  begins  with  a  description  of  the  scattering  matrix  and  a  summary  of 
the  CPM  in  Section  2.  This  section  is  designed  for  those  who  have  studied  basic  quantum 
mechanics  but  have  had  no  exposure  to  the  CPM.  It  discusses  a  derivation  of  the 
scattering  matrix,  the  Moller  operators,  the  split-operator,  how  to  formulate  the  reactant 
and  product  wave  packets,  and  the  approach  to  calculate  the  S-matrix  from  the  correlation 
functions.  I  attempt  to  make  the  discussion  clearer  by  explaining  how  these  topics  are 
applied  to  a  square  well  scattering  problem.  In  section  3,  the  requirements  of  the  previous 
algorithm  developed  in  [1-3]  are  discussed  as  well  as  explanations  of  why  the  algorithm 
fails  at  low  energy. 

The  new  algorithm  developed  in  this  thesis  is  then  presented  with  a  focus  on 
discussing  how  it  will  result  in  significant  improvements  in  the  low  energy  regime. 
Section  4  then  describes  the  parameters  of  a  typical  CPM  simulation  as  well  as  the 
implications  of  changing  them.  In  section  5  the  results  of  low  energy  scattering  are 
presented  for  a  square  well.  Some  of  the  difficulties  encountered  are  discussed  as  well  as 
the  convergence  rate  of  the  error  associated  with  a  discrete  square  well  based  on  the 
coordinate  grid  spacing.  Section  6  discusses  the  first  of  two  error  analyses;  the 
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methodology  behind  the  error  analysis  is  presented  along  with  the  results.  Seetions  7  and 
8  diseuss  the  trapezoidal  well  as  well  as  the  seeond  error  analysis  showing  that  the  S- 
matrix  elements  produeed  by  the  new  algorithm  eonverge  to  the  eorreet  solution  in  the 
low  energy  regime.  Finally,  a  summary  of  the  thesis  is  presented  as  well  as 
reeommendations  for  future  researeh. 
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II  Time  Dependent  Scattering  Theory 

2,1  The  S-Matrix 

As  noted  earlier  the  seattering  matrix  is  related  to  the  transmission  and  refleetion 
ooeffieients  introdueed  in  a  basie  quantum  meehanies  eourse.  To  demonstrate  the 
relationship  I  will  adopt  a  problem  from^^n  Introduction  to  Quantum  Mechanics  by 
David  J.  Griffiths  [4],  Consider  an  arbitrary,  loealized  potential 

0  x<-a 

V{x)  =  <V{x)  -a<x<  a,  (1) 

0  x>  a 


y/{x)  =  Fd'°‘  +Ge-''°‘ .  (3) 


Region  I  Region  II  Region  III 

Figure  1 .  Scattering  from  an  arbitrary  localized  potential  (  V (x)  =  0  except  in  Region  II) 
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In  Region  II,  \j/  cannot  be  determined  until  the  potential  is  specified,  but  because  the 
Schrodinger  equation  is  a  linear,  second-order  differential  equation,  the  general  solution 
has  to  be  of  the  form 

W{x)  =  Cf{x)  +  Dg{x),  (4) 

where  / (x)  and  g(x)  are  two  linearly  independent  particular  solutions  [5].  There  will  be 
four  boundary  conditions: 

Wx  i-a)  =  Wi  (-«) 

\j/^{a)=\j/^{a) 
y/\{a)  =y/\{a) 

where  the  subscript  indices  1,  2,  and  3  stand  for  the  three  regions.  Two  of  these  boundary 
conditions  can  be  used  to  eliminate  C  andZ) ,  and  the  other  two  can  be  “solved”  for 
B  and  F  in  terms  of  A  and  G  : 


B  =  S,,^_,A  +  S^,^_,G, 

(6) 

and 

(V) 

The  four  coefficients 5*+^  depend  on  k  (and  hence £■)  and  constitute  the  2x  2  S  -matrix. 


The  S-matrix  tells  you  the  outgoing  amplitudes  (B  and  F)  in  terms  of  the  incoming 
amplitudes  ( A  and  G  ): 


(8) 


Considering  only  scattering  from  the  left  ( G  =  0  ); 
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|2 


(9) 


Ri  = 


\B- 

A\ 


=  S_ 


+k,-k\  ’ 


G=0 


F 

u 


=  s_ 


G=0 


or  from  the  right  (A  =  0), 


R..  = 


IfP  I  ,2 

B^  1  i 

II  —  C 

1  ,2  ’ 

1  A=0 

A=0 

(10) 


where  the  subseripts  /  and  r  denote  seattering  from  left  and  right  respeetively.  Any 
realistie  seattering  experiment  will  eonsist  of  sending  a  beam  of  partieles  into  the 
interaetion  region.  The  beam  will  have  some  spread  in  energy  (henee  momentum)  and 
thus  eannot  be  represented  as  a  plane  wave.  To  handle  the  spread,  Gaussian  wave 
paekets  are  used  to  represent  the  (inbound)  and  (outbound)  states.  In  the  more  realistie 
pieture,  equation  (8)  will  take  on  a  more  eomplieated  form  eonsisting  of  a  “bloek” 
diagonal  S-matrix, 


f  t  ^ 

^  \ 

\ 

r  t  ^ 

F{+k) 

9  9 

^+k,+k  ^-k,+k 

A{+k) 

B{-k) 

G{-k) 

F{+k') 

9  9 

^+k\+k’  ^-k\+k' 

A{+k') 

B{-k') 

Gi-k') 

[  ^  J 

1 

\ 

) 

i 

\  ^  J 

where  the  arrows  indieate  that  k  e  |-oo,oo|  and  the  expansion  eoeffieients  A,B,G,  and  F 


represent  the  Fourier  transform  of  the  reaetants  and  produets  [6].  For  simple  potentials, 
sueh  as  a  square  well,  solving  the  Sehrodinger  equation  in  the  interaetion  region  is 
straightforward  and  the  S-matrix  elements  ean  be  obtained  analytieally.  More  interesting 
and  realistie  potentials  however,  ean  beeome  diffieult  or  impossible  to  solve  analytieally 
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and  require  numerieal  teehniques.  Both  time  dependent  and  time-independent  approaehes 
have  been  derived  to  eompute  S-matrix  elements  numerieally.  The  Channel  Paeket 
Method,  a  time-dependent  approaeh,  is  used  to  eompute  the  S-matrix  elements  in  this 
thesis. 

2.2  The  Channel  Packet  Method 

As  noted  earlier  the  CPM  was  developed  by  Dr.  David  Weeks  and  Dr.  David 
Tannor  [1-3].  The  CPM  is  a  proeess  through  whieh  S-Matrix  elements  are  ealeulated 
numerieally  by  applying  the  elassieal  Moller  operators  (13)  to  propagate  the  reaetant 
states  I  'P,^ )  and  the  produet  states  |  'P^^, ) .  The  “seattering”  equation  relates  the  ineoming 
state  with  the  outgoing  state: 

k„,>=s|r,.>,  (12) 


where  S  denotes  the  “seattering”  operator.  The  CPM  formulates  the  S-matrix  elements 
that  make  up  S  through  a  time-dependent  eorrelation  funetion  between  what  are  known 
as  Moller  states.  The  Moller  states  are  formed  by  applying  the  Moller  operators  to  the 
reaetant  and  produet  states. 

2,2.1  Moller  Operators 

The  Moller  operators  are  defined  as: 

=  lim  [Qxp{iHt  /  h)  exp{-iHj  /h)],  (13) 

?->-+0D 


where  Ho  and  // represent  the  asymptotie  and  full  Hamiltonians: 

2 


P 

2m  ’ 


(14) 
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and 


(15) 


//  =  ^  +  F(x). 
2m 


The  relationship  between  the  full  and  asymptotie  Hamiltonians  is  defined  as: 


(16) 


The  operators  are  related  to  the  seattering  operator  by  [7] 

S  =  0}  Q  ,. 


(17) 


The  relationship  defined  by  equation  16  requires  that  the  potential  go  to  zero  as  x  ^  ±oo  . 
The  Moller  operators  aet  on  the  reaetant  and  product  states  to  form  the  Moller  states  in 
the  following  way: 


y) 

Win)^ 

(18) 

and  1 

¥-) 

=  Q_ 

Won)- 

(19) 

The  plus  and  minus  superscripts  denote  the  “reactant”  and  “product”  Moller  states 
respectively.  Numerically  applying  the  Moller  operators  results  in  the  reactant  state 
I  ) ,  being  propagated  to  t  =  -go  under  the  asymptotic  Hamiltonian  and  then  propagated 

to  t  =  0  under  the  full  Hamiltonian.  Conversely,  the  product  state  ,  is  propagated  to 

t  =  GO  under  the  asymptotic  Hamiltonian  and  then  propagated  to  t  =  0  under  the  full 
Hamiltonian.  Numerically  speaking,  the  states  cannot  be  propagated  to  t  =  ±go  ,  however, 
the  asymptotic  limits  are  achieved  when  the  states  have  propagated  completely  outside 
the  interaction  region  [9].  To  better  understand  the  application  of  the  Moller  operators  a 
step-by-step  discussion  is  provided  below  for  the  specific  case  of  a  square  well  potential. 
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Figures  2  and  3  show  the  eoordinate  and  momentum  representations  of  a  typieal  CPM 
setup.  Gaussian  wave  paekets  represent  the  produet  and  reaetant  states.  The  eoordinate 
representation  of  an  initial  state  is  defined  as: 

0)  =  1/2,0  v/4  (x-xj]  exp  [-(x  -  X /  4cr/  ]  ,  (20) 

cr^  (2;r) 

where  cr^  deseribes  the  initial  spread  in  the  position  of  the  partiele  and  serves  to  define 
the  eenter  of  the  wave  paeket.  The  eorresponding  probability  density  takes  the  form  of  a 
Gaussian  funetion  and  is  deseribed  as: 

T-f  =— t^exp[-(x-x„)V2<T/].  (21) 

a^yJZTT 


Figure  2.  The  initial  product  and  reactant  Gaussian  wave 

packets  I ^  and  )  plotted  in  the  coordinate  representation. 

The  square  well  potential  depth  is  scaled  by  a  factor  of  10.  The 
probability  density  is  shown. 
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Figure  3.  The  initial  product  and  reactant  Gaussian  wave 
packets  I ^  and  1^^^^  ^  plotted  in  the  momentum  representation, 
where  only  ’\-k  is  considered.  The  probability  density  is  shown. 


The  complex  modulation  exp[zA:^(x-x_,)]m  equation  (20)  serves  to  define  the  average 
momentum 

(p)  =  trK.  (22) 

The  physical  interpretation  of  the  state  defined  by  equation  (20)  is  then  a  particle 
localized  within  a  spread  of  <j^  about  moving  with  an  average  momentum,  .  The 

expansion  coefficients  of  the  initial  states  are  then  computed  using  the  Fourier  transform 
of  y/{x,Q) : 


ri{k)  =  -^\  'F(x,0)exp( 


-ikx)dx  , 


with  the  corresponding  probability  density 

2cr 


dikf  =^=Qx^[-2cT^\k^  -kf^ 
■\J27r 


(23) 


(24) 
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The  momentum  probability  density  is  thus  also  a  Gaussian  eentered  about  with  a 


spread  cr^  =  (2cr^)  ' .  The  produet  of  the  uneertainties  has  its  minimum  value  in  the  initial 
Gaussian  wave  paeket: 


h 


Gauss 


2  .  (25) 

As  diseussed  earlier  the  states  must  be  propagated  under  the  asymptotie  Hamiltonian 
during  the  eonstruetion  of  the  Moller  states.  Sinee  V{x)  =  0  in  the  asymptotic  regions,  the 
states  propagate  as  free  particles  during  the  first  half  of  the  Moller  operators.  The 
blessing  in  this  is  that  one  can  propagate  to  the  asymptotic  channels  analytically  using  the 
time-dependent,  free-particle,  Gaussian  state  [8]; 


exp 


.  T 


V  2f^.  y 


exp 


-{ir  /  4tcr^)(x  -x^-  hkj  /  mY 
l  +  it  /  T 


cr. 


1/2 


(26) 


where, 

laYm 

T= - ^ - 

n 


(27) 


However,  in  doing  so  there  is  a  distortion  of  the  Gaussian  states  in  time.  That  is,  the 
initial  uncertainty  in  position,  <jYK)  ’  becomes  larger  in  time  and  the  wave  function 
spreads  as 


2  +  2 


1  +  - 


t  n 


,1/2 


[8]. 


(28) 
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The  spreading  can  be  understood  by  realizing  that  the  initial  Gaussian  state  is  a 
superposition  of  plane  waves,  each  having  different  momentum.  Thus,  as  time  progresses 
the  initial  uncertainty  in  position  is  compounded  by  the  initial  uncertainty  in  the 
momentum. 

The  coordinate  representation  (Figure  2)  shows  two  Gaussian  wave  packets 
centered  in  the  middle  of  a  square  well  potential.  The  reactant  and  product  states  are 
constructed  to  be  the  same  in  this  case.  The  results  of  the  application  of  the  first  part  of 
the  Moller  operators  are  shown  in  Figure  4.  The  reactant  wave  packet  has  propagated 
backwards  in  time  300000  atomic  units  (a.u.)  (to  the  left  since  it  has  positive 
momentum),  and  the  product  has  propagated  forward  in  time  300000  a.u.  (to  the  right 
since  it  too  has  positive  momentum).  The  propagations  are  performed  analytically  using 
equation  (26).  These  states  are  defined  as  the  intermediate  Moller  states 

(29) 


□.3 

□  2S 

□.2 

□  IS 

□.1 

□  US 

V(x) 

-□US 

■20  -15  -ID  -S  □  5  ID  15  21 

X  {atomic  units} 

Figure  4.  The  probability  density  of  the  intermediate  Moller  states, 
and  (after  t  -  30000  a. u. )  plotted  in  the  coordinate 

representation.  The  square  well  potential  is  also  plotted. 


=  lim  exp(-///^t/  fi) 

/  /— »+oo 


'F. 
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Both  propagations  have  lasted  long  enough  to  ensure  the  wave  paekets  are  eompletely 
outside  the  square  well.  It  should  be  noted  that,  by  eonstruetion,  these  states  eontained 
only  positive  momentum  (Figure  3).  Therefore,  it  is  a  simple  matter  to  analytieally 
propagate  the  states  eompletely  outside  the  interaetion  region. 

The  final  step  in  the  eonstruetion  of  the  Moller  states  is  to  apply  the  seeond  half  of 
the  Moller  operators  to  these  intermediate  states.  This  applieation  is  mathematieally 
defined  as: 


y/^')=  lim  exp(z//t//z) 
/  ?— >+00 


where  H  eontains  the  square  well  potential  (atomie  units), 

f-0.01  -2<x<2 


V{x)  =  - 


0 


elsewhere 


(30) 


(31) 


The  applieation  of  equation  (30)  is  done  numerieally  using  the  split  operator.  The  result 
of  this  applieation  is  shown  in  Figure  5.  The  states  have  propagated  baek  to  t  =  0 .  Sinee 
the  produet  and  reaetant  eontain  only  positive  momentum,  only  S^^.  ( the  probability  of 

eoming  in  with  positive  momentum  and  going  out  with  positive  momentum)  may  be 
eomputed.  The  Moller  operators  are  speeifie  to  seattering  theory,  the  ehoiee  of 
propagation  method,  however,  is  not.  In  this  thesis  the  method  of  propagation  is  the  split- 
operator  method. 

2.2,2  The  Split  Operator  Method 

The  Split  Operator  is  used  in  the  CPM  to  propagate  the  initial  produet  and 
reaetant  wave  funetions  mentioned  in  the  previous  seetion. 
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0J6 


Figure  5.  The  M0ller  states,  |t+)  and  |t-)  plotted  in  the  coordinate 
representation.  The  square  well  potential  has  been  scaled  by  a  factor  of 
10.  The  probability  density  is  shown. 

The  split  operator  approximates  the  time  evolution  operator, 


t/  =  exp  Z(T-\-V) 


(32) 


where, 


3  •  ^ 

A  =  -i  — 

n 


(33) 


as 


U  »  exp 

( 

A- 

0 

exp| 

jexp 

(  f  ^ 
A- 

0 

1 

1 

(34) 


The  splitting  eontains  an  error  sinee  the  eommutator,  \y{x),T~\  ^  0 ;  the  error  term  ean  be 
easily  obtained  from  a  Taylor  expansion  of  the  exponentials  (Appendix  A).  Terms 
appearing  in  the  split  operator  are; 
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(35) 


and 


T  = 


2m 


V{Xj), 


(36) 


(37) 

(38) 


The  index  j  is  defined  over  the  set  described  in  equation  (38)  specifically  for  the 
application  of  using  the  Fast  Fourier  Transform  (FFTs)  where  is  the  number  of  grid 
points.  Propagation  is  achieved  by  applying  the  approximate  time  evolution  operator 
listed  in  equation  (34)  to  the  product  and  reactant  wave  packets.  Using  Fast  Fourier 
Transforms  one  can  rapidly  switch  between  the  position  and  momentum  representations 
of  the  reactant  and  product  wave  packets.  The  time  evolution  operator  U  is  split  into 
position  and  momentum  dependent  parts.  The  advantage  of  using  the  approximate,  split 
operator  is  that  the  position  and  momentum  parts  of  the  operator  become  simple 
multiplicative  factors  in  the  coordinate  and  momentum  representations.  The  latter  point 
saves  computational  effort.  To  explain  this  point,  consider  applying  the  time  evolution 
operator  to  a  Gaussian  wave  packet  |'F(x)) .  In  order  to  characterize  the  evolving  wave 
packet  one  would  use; 


‘F(x.,t.^i)  =  exp 


2m  dx^ 


+  F(x,) 


Wix„0 


(39) 
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Each  time  step  requires  the  use  of  a  finite  differenee  method  resulting  in  a  substantial 
increase  in  the  total  eomputational  effort  needed  to  fully  propagate  the  wave  packets. 
Conversely,  with  the  use  of  the  split  operator  method,  a  EFT  is  applied  to  ‘E(x)  to  form 
‘E(A:)  .  The  differential  representation  of  the  momentum  operator  is  avoided,  sinee  in  the 
momentum  representation,  the  operator  is  diagonal.  The  process  of  applying  FFT’s  and 
IFFT’s  (FFT"')  to  switeh  baek  and  forth  between  the  eoordinate  and  momentum 
representation  eontinues  with  a  single  time  step  consisting  of 

r  r  r  r  1T 

^^{x,t^^^)  =  IFFT  exp  X-  FFT  exp(;LF)/FFr  exp  X-  FFT[^>{x,t,)]  .(40) 

v^/L  Lv^y  JJ_ 

The  proeess  is  continued  over  and  over  again  until  the  desired  propagation  time  is 
achieved.  The  application  of  the  split  operator  is  used  to  propagate  the  reactant  and 
product  wave  paekets  in  all  simulations  performed  in  this  thesis.  The  simulations 
performed  do  not  aecount  for  the  error  assoeiated  with  the  split  operator.  However,  At  is 
made  suffieiently  small  (1  a.u.)  in  an  effort  to  make  the  error  negligible. 

2,2,3  Computing  the  S-Matrix  Elements 

After  the  Gaussian  produet  and  reaetant  wave  paekets  are  formulated  and  the 
Moller  states  construeted,  the  S-matrix  elements  are  computed.  A  convenient  starting 
point'  for  diseussing  how  we  attempt  to  ealeulate  the  S-matrix  elements  starts  with  the 
Fourier  Transform  (FT)  of  the  correlation  function  between  the  reaetant  and  produet 
Moller  states; 

*  For  the  derivation  of  this  equation  see  D.E.  Weeks  and  D.J.  Tannor,  Chem.  Phys.  Letters  207  (1993)  301. 
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I-  1  roo  ATTffl 

FT[C{t)]  =  dt  exp(/£0C(0  =  ^  ^  (+^) 

+  r]-{+k)r]\+k)  ,  (41) 

+  ri-{-k)ri\-k)  S_,  _, 

+  7J-{+k)7J^{-k) 

where, 

C(0  =  (t-  |exp(-///0|'I'^)  .  (42) 

and  the  plus  and  minus  superseripts  in  the  7“  stand  for  reaetant  and  produet  respeetfully. 
The  expansion  eoeffieients  7j~(±k)  listed  in  equation  (41)  are  the  same  as  those  defined 
previously  in  equation  (23).  The  plus  and  minus  designations  in  the  S-matrix  elements 
stand  for  the  four  possible  eombinations  in  a  seattering  model.  For  example,  , 

would  deseribe  the  S-matrix  element  for  an  ineoming  state  (reaetant)  having  positive 
momentum  and  an  outgoing  state  (produet)  leaving  with  negative  momentum.  The 
absolute  value  of  squared  would  then  give  the  probability  for  eoming  in  from  the 

left  with  positive  momentum  and  leaving  to  the  left  with  negative  momentum.  The  time- 
dependent  eorrelation  funetion  of  equation  (42)  is  numerieally  evaluated  using  the  split- 
operator  by  propagating  using  the  full  Hamiltonian  until  it  is  eompletely  outside  the 

interaetion  region.  At  eaeh  time  step  the  sealar  produet  given  in  equation  (42)  is 
eomputed.  It  is  important  to  note  that  in  order  to  eompletely  form  the  eorrelation 
funetion  one  must  propagate  |  ^  both  forward  and  baekward  in  time  until  |  ^  is 

eompletely  outside  the  interaetion  region. 
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Ill  Application  of  the  CPM 


3.1  Previous  Technique 

Past  approaches,  such  as  those  described  in  [1-3],  have  been  to  simply  invert 
equation  (41)  to  establish  the  S-  matrix  elements, 

^  i»iWi-nc(o] 

2nm  rj  {±k)ri^{±k) 

In  order  to  write  equation  (43)  specifie  requirements  must  be  imposed  on  the  product  and 
reactant  states.  That  is,  the  produet  and  reactant  Gaussian  wave  packets  must  contain  only 
positive  or  negative  momentum  ranges.  For  example,  if  one  wishes  to  eompute  the  S- 
matrix  element  ,  the  reactant  and  product  wave  packets  must  contain  only  positive 

momentum  ranges  (Figure  3).  Equation  (41)  then  simplifies  to 


\fi\\k\FT[C{t) 

2nm  ri  {+k)ri^  {+k)  ’ 


(44) 


since  the  r]  {-k)  and  r]^{-k)  terms  all  go  to  zero.  Similarly,  the  other  S-matrix  elements 
may  be  computed  under  the  guidelines  given  in  table  1.  Weeks  and  Tannor  developed 
this  method  to  formulate  S-matrix  elements  with  excellent  results.  The  limitation  in  using 
there  approach  however,  was  that  the  product  and  reactants  could  only  contain  either 
positive  or  negative  momentum.  As  a  result  they  could  not  accurately  ealculate  S-matrix 
elements  near  zero  energy. 

The  explanation  why  previous  attempts  failed  to  produce  good  results  at  low 
energies  is  that  there  was  not  enough  information  eontained  in  the  initial  wave  packets. 
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Table  1.  S-matrix  selection  table. 


To  Compute 

Use 

r]{k)  =  0 

A:<0 

k>0 

k<0 

A:>0 

rj^ik)  =  0 

A:<0 

A:<0 

k>0 

A:>0 

And  divide 

FT[C{t)]hy 

ri{+k)ri^  {+k) 

ri{+k)ri\-k) 

ri{-k)ri\+k) 

ri{-k)ri\-k) 

The  requirement  that  the  reactant  and  product  states  contain  only  positive  or  negative 
momentum  leads  to  a  loss  of  information  at  low  energies.  Consider  the  product  and 
reactant  states  shown  in  Figure  3.  The  wave  packets  are  centered  about  =  6 .  Notice 

that  the  probability  density  near  the  low  momentum  ranges  is  effectively  zero.  Of  course, 
the  theoretical  wave  packet  extends  to  plus  and  minus  infinity,  and  thus  would  be  fully 
represented  across  the  entire  momentum  range.  However,  computers  are  limited  in 
storing  very  small  numbers  and  thus  the  amplitude  in  the  low  momentum  range  is 
interpreted  as  zero;  leading  to  a  loss  of  information  for  approximately  k  <  1 .  It  should  be 
noted  as  well  that  another  loss  of  information  is  incurred  for  approximately  A:  >  1 1 .  The 
obvious  solution  is  to  center  the  wave  packets  closer  to  zero. 

Consider  the  wave  packet  shown  in  Figure  6.  Notice  that  there  is  a  substantial 
probability  density  near  low  momentum.  However,  as  the  center  of  the  wave  packets  are 
moved  into  the  low  momentum  regime  the  wave  packets  overlap  zero  and  exhibit 
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substantial  amplitude  in  the  negative  momentum  regime.  The  requirement  that  the 
produet  and  reaetant  wave  packets  contain  only  positive  or  negative  momentum  ranges  is 
violated. 


k  {atomic  units} 

Figure  6.  Initial  product  and  reactant  Gaussian  wave  packets 
plotted  in  the  momentum  representation  showing  overlap.  The 
probability  density  is  shown. 

Previously,  Weeks  and  Tannor  simply  ignored  the  overlapping  amplitude  and  proceeded 
with  equation  (43)  under  the  formulations  made  in  Table  1.  In  an  effort  to  reduce  the 
error  associated  with  neglecting  the  overlapping  momentum  amplitude  they  centered  the 
wave  packet  such  that  there  was  as  little  overlap  as  possible.  However,  this  becomes 
problematic  since  one  must  carefully  construct  the  wave  packets  such  that  the  tails  go  to 
zero  at  A:  =  0  .  In  doing  so  however,  one  is  always  fighting  the  battle  of  having  the  tails  go 
to  zero  while  still  needing  information  to  be  numerically  stored  at  low  energies. 

A  new  problem  is  introduced  if  the  reactant  or  product  state  has  both  positive  and 
negative  momentum.  Imagine  applying  the  Moller  operator  to  a  reactant  state  centered 
in  the  middle  of  a  potential  well  like  that  of  Figure  1  with  a  momentum  representation 
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like  that  of  Figure  6.  As  the  wave  packet  is  propagated  backwards  in  time  under  the 
asymptotic  Hamiltonian,  most  of  the  wave  packet  will  move  to  the  left,  however  some 
would  move  to  the  right  since  it  contains  some  negative  momentum.  The  time  required  to 
fully  propagate  the  reactant  and  product  wave  packets  completely  outside  the  interaction 
region  increases  due  to  movement  in  both  directions. 

The  subsequent  increase  in  propagation  time  will  require  a  larger  grid  since, 
according  to  according  to  equation  (28),  the  wave  packet  will  spread  as  a  function  of 
time.  Since  Weeks  and  Tannor  used  wave  packets  with  all  positive  or  negative 
momentum  ranges  this  problem  was  avoided;  it  is  fairly  easy  to  propagate  the  wave 
packets  outside  the  interaction  region  when  they  only  contain  positive  or  negative 
momentum.  In  order  to  apply  the  CPM  to  wave  packets  in  the  low  energy  regime  one 
must  reformulate  the  previous  application  and  loosen  the  restriction  that  the  product  and 
reactant  wave  packets  contain  only  positive  or  negative  momentum. 

3,2  New  Technique 

In  order  to  establish  enough  information  in  the  low  energy  regime,  one  formulates 
the  product  and  reactant  wave  packets  to  have  both  positive  and  negative  momentum.  In 
doing  so  however,  the  requirements  allowing  the  use  of  equation  (43)  are  violated. 
Therefore,  a  new  approach  must  be  applied  to  the  CPM  to  handle  the  calculation  of  S- 
matrix  elements  at  low  energy.  The  algorithm  developed  in  this  thesis  extends  equation 
(41)  into  a  matrix  formulation  with  four  product  and  reactant  wave  packets.  The 
mathematical  matrix  equivalent  is 
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where  the  rjj  -matrix  is  given  by, 

'  ri-{-k\ri\+k\  rj- {+k\ri\+k\  rj- {-k\ri\-k\  rj- {+k\ri\-k\ 

ri  {-k)^ri\+k)^  ri  {+k)^ri\+k)^  ri  {-k)^ri\-k)^  ri  {+k)^ri\-k)^ 
ij!  =  .  (46) 

V  (-^)3^^(+^)3  V  V  (-^)3?7^(-^)3  V 

ri-{+k)^ri\+k)^  rj- {-k) ,ri\-k) ,  ri-{+k)^ri\-k),^ 

Using  a  four-by-four  matrix  approach  one  is  able  to  handle  product  and  reactant  wave 
packets  containing  any  combination  of  positive  and  negative  momentum.  The  entire  S- 
matrix  is  then  computed  by  inverting  equation  (45); 

V...  -[C(0,] 

Ittm  Tr[C(03] 

U-.-J  ^Tr[C(04], 

In  the  limit  that  the  product  and  reactant  wave  packets  are  constructed  to  have  only 
positive  or  negative  momentum,  the  rji  -matrix  defined  by  equation  (46)  becomes 
diagonal,  and  reduces  to  the  four  independent  S-matrix  equations  of  (43). 

The  requirement  of  equation  (47)  is  that  the  inverse  of  rji  must  exist.  Therefore,  a 
necessary  condition  is  that  the  rji  -matrix  must  have  a  non-zero  determinant.  Numerically 
representing  a  Gaussian  wave  packet  will  always  lead  to  expansion  coefficients  that  are 
interpreted  as  zero  near  the  tail  end  of  the  wave  packet.  In  turn,  the  rjj  -matrix  will  be 
fdled  with  zeros  and  therefore  become  singular.  However,  as  in  previous  approaches,  this 
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is  simply  interpreted  as  an  inability  to  extraet  information  from  the  tail  ends  of  the 
produet  and  reaetant  wave  paekets.  Therefore,  ;^need  only  to  be  non-singular  in  the 
energy  range  of  interest. 

Due  to  the  nature  of  a  four-by  four  matrix  approaeh  one  must  formulate  four 
separate  product  and  reactant  states.  One  could  simply  establish  four  wave  packets 
centered  at  four  distinct  momentum  values.  In  doing  so,  one  would  ensure  the  t!  -matrix 
is  non-singular  in  the  energy  range  of  interest.  In  turn,  four  propagations  would  be  needed 
to  formulate  the  four  transforms  listed  in  equation  (47).  However,  if  one  is  careful  in 
their  construction  of  the  product  and  reactant  states,  only  two  propagations  are  required. 
To  explain,  consider  formulating  the  product  and  reactant  wave  packets  in  the  following 
way.  Establish  two  Gaussian  wave  packets,  one  centered  about  and  the  other 

about .  Reactant  1  is  made  to  be  the  same  as  product  1  and  product  2  the  same  as 

reactant  1 .  The  choice  of  where  to  put  the  wave  packets  in  the  coordinate  representation 
is  arbitrary.  The  setup  is  shown  in  Figures  7  and  8. 

Given  the  initial  setup  described  above  one  then  numerically  computes  the 
corresponding  Moller  states  for  each  product  and  reactant.  Four  correlation  functions  are 
then  computed  to  establish  the  right  hand  side  of  equation  (47).  Given  the  setup 
described  above,  the  four  correlation  functions  would  be 
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(48) 


C,(0  =  ('F-Jexp(-///0|'F\) 
Q(0  =  ('F-Jexp(-///0|'F"_) 
C3(0  =  ('F-Jexp(-///0|'F^J 
Q(0  =  ('F-Jexp(-///0|'F^) 


Figure  7.  Initial  product  and  reactant  Gaussian  wave 


packets , 


plotted  in  the  coordinate  representation.  The 


square  well  potential  has  been  scaled  by  a  factor  of  10.  The 
probability  density  is  shown. 


k  {atomic  units} 

Figure  8.  Initial  product  and  reactant  Gaussian  wave  packets  \  and 
plotted  in  the  momentum  representation.  The  probability  density  is  shown. 
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Once  the  eorrelation  funetions  are  eomputed  a  simple  Fourier  transform  is  performed  and 
the  S-matrix  is  determined  using  equation  (47). 

3,2.1  Summary  of  Steps  for  New  Approach 
A  summary  of  the  steps  used  to  perform  a  simulation  is  listed  below. 

1 .  Make  produet  1  and  reaetant  1  equal  to  one  another  and  eontain  nearly  all  positive 
momentum. 

2.  Make  produet  2  and  reaetant  2  equal  to  one  another  and  eontain  nearly  all 
negative  momentum. 

3.  Using  equation  (23)  fill  in  the  rjj  -matrix  based  on  the  produet  and  reaetant  wave 
paekets  from  steps  1  and  2  as  shown  below.  The  subseripts  1  and  2  represent 
reaetant  and  produet  one  and  two. 

'  ri  {-k\ri\+k\  ri  {+k\ri\+k\  Ti  {-k\Tj\-k\  ri  {+k\ri\-k\  ^ 
ri  {-k\ri\+k\  ri  {+k\ri\+k\  ri  {-k\ri\-k\  ri  {+k\ri\-k\ 

^  ri  {-k\7i\+k\  ri-{+k\ri\+k)^  rj- {-k\7i\-k)^  7i  {-vk\ri\-k\ 
{-k\ri\+k\  rj- {+k\ri\+k\  rj- {-k\ri\-k\  ri  {+k\ri\-k\ ^ 

4.  Apply  the  split  operator  to  numerieally  eompute  the  Moller  states  for  eaeh  produet 
and  reaetant  under  the  operations  defined  in  equations  (18)  and  (19). 

5.  Use  the  split  operator  to  propagate  the  reaetant  Moller  states  both  forward  and 
baekward  in  time  under  the  full  Hamiltonian  while  eomputing  the  eorrelation 
funetions  defined  by  equation  (48)  at  eaeh  time  step. 
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6.  Compute  the  Fourier  transforms  to  eonvert  the  time  dependent  eorrelation 
funetions  into  to  the  energy  representations  listed  on  the  right  hand  side  of 
equation  (47). 

7.  Finally  use  equation  (47)  to  eompute  the  S-matrix  elements. 

Using  the  algorithm  deseribed  above  one  ean  eompletely  eharaeterize  wave 
paekets  eontaining  both  negative  and  positive  momentum  ranges.  By  keeping  the  full  rjj  - 
matrix  information  ean  be  stored  in  the  low  energy  regime;  thus  improving  the 
ealeulations  of  the  S-matrix  elements  aeross  the  entire  energy  speetrum.  Furthermore,  the 
method  used  is  more  eomplete  than  previous  methods  sinee  it  does  not  require  reaetant 
and  produet  states  having  all  positive  or  negative  momentum. 

3.2,2  An  Interesting  Implication  of  the  New  Approach 

It  should  be  noted  that  aside  from  the  full  matrix  approaeh  to  the  CPM  there  is  an 
interesting  differenee  between  equations  (47)  and  (43).  Equation  (47)  implies  that  when 
reaetant  and  produet  states  eontain  both  positive  and  negative  momentum,  one  eannot 
know  one  S-matrix  element  without  simultaneously  knowing  the  other  three.  That  is, 
there  is  an  inter  relationship  among  the  S-matrix  elements.  Equation  (43)  requires  no 
eonnection  whatsoever  between  the  different  S-matrix  elements  exeept  that  the  sum  of 
the  transmission  and  refleetion  eoeffieients  be  one.  Eurthermore,  in  the  limit  that  the 
produet  and  reaetant  wave  paekets  eontain  only  positive  or  negative  momentum  ranges, 
the  -matrix  beeomes  diagonal  and  the  S-matrix  elements  beeome  “detangled.” 
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IV  Choice  of  Parameters 


There  are  several  parameters  relevant  to  the  CPM.  These  parameters  inelude: 
the  grid  spaeings  Ax,  AA:,  At,  and  AE ,  the  eentering  x„  and  of  the  Gaussian  produet  and 

reaetant  wave  paekets,  the  maximum  and  minimum  values  of  the  grid  sizes  , 

max  max 

and  ,  the  width  of  the  eoordinate  and  momentum  wave  paekets  cr^  and  ,as  well  as 

max 

the  number  of  grid  points  N  .  Many  of  the  parameters  mentioned  are  determined  using 
knowledge  about  the  setup  of  the  seattering  problem  itself  The  parameter  x^  is 

somewhat  arbitrary,  and  in  the  simulations  performed  in  this  thesis,  was  ehosen  to  be  the 
origin.  The  eentering  of  the  momentum  wave  paeket  however,  has  a  mueh  more  profound 
effeet.  The  parameters  and  cr^.  will  speeify  the  range  of  energy  for  whieh  S-matrix 

elements  can  be  computed.  Due  to  the  nature  of  this  thesis,  this  parameter  is  chosen  such 
that  the  momentum  representation  contains  both  positive  and  negative  momentum.  The 
product  and  reactant  Gaussian  states  described  in  equation  (20)  require  a  specification  for 
the  initial  width  of  the  wave  packet.  One  thing  to  consider  in  choosing  the  parameter  cr^ 

is  that  the  smaller  it  is  the  wider  the  wave  packet  becomes  in  the  momentum 
representation.  A  wide  momentum  wave  packet  corresponds  to  the  ability  to  calculate  S- 
matrix  elements  over  a  wide  range  of  energies  during  a  given  simulation.  However,  the 
wider  the  wave  packet  becomes  the  bigger  the  grid  size  needed  to  fully  represent  it. 

Perhaps  the  most  important  computational  parameter  is  N  .  Since  FFT’s  are  used  to 
propagate  the  wave  packets  this  parameter  must  satisfy 

A  =  2”,ne{l,2,3...}.  (49) 
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As  the  number  of  grid  points  is  inereased  more  computational  effort  is  needed  to  perform 
each  FFT.  Once  N  is  specified  the  grid  spacing  Ax  is  computed  using: 

AX=  "^max 

N 

The  parameters  x^^^  and  are  determined  by  considering  the  grid  size  needed  to 
accurately  represent  all  states  during  the  simulation.  The  intermediate  Moller  states,  like 
those  shown  in  Figure  4,  define  the  states  for  which  the  wave  packets  are  as  wide  and  as 
far  out  as  any  state  during  a  simulation.  Hence,  they  can  be  used  as  limiting  cases  in 
determining  the  required  coordinate  grid  size.  A  further  consideration,  specific  to  a 
square  well,  is  to  ensure  that  Ax  is  made  small  enough  such  that  the  discontinuities 
associated  with  the  potential  are  sufficiently  represented.  To  illustrate  this  point  consider 
the  two  potentials  shown  in  Figure  9.  These  diagrams  are  meant  to  show  the  effect  of 
decreasing  the  coordinate  grid  spacing.  A  large  Ax  results  in  a  more  trapezoidal  figure. 
Thus,  in  an  effort  to  represent  the  square  well  one  must  ensure  the  coordinate  grid 
spacing  is  sufficiently  small. 


Figure  9.  (a)  Approximate  square  well  potential  for  large  Ax  .  (b)  Approximate  potential 
showing  the  convergence  to  an  accurate  square  well  as  Ax  is  decreased. 
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Since  FFTs  are  used  to  switch  back  and  forth  between  the  eoordinate  and  momentum 
representations,  the  momentum  grid  spaeing  is  determined  by  a  eharaeteristie  of  the  FFT. 
The  FFT  requires  a  speeial  relationship  between  the  funetions  transformed  from  one 
representation  to  the  other.  The  relationship  is 


= 


In 

NAx  ' 


(51) 


The  parameters  and  naturally  follow  as 


k 

max 

min 


(52) 


where  kj  =  jAk,  2  ^  y  ^  yj  •  (53) 

The  size  of  the  momentum  grid  must  be  able  to  sustain  the  momentum  wave  paekets 
during  the  entire  simulation.  A  good  estimate  for  ensuring  that  the  momentum 
representations  will  be  defined  throughout  the  entire  simulation  is  to  perform  a  quiek 
estimate  of  the  maximum  value  of  momentum  aehieved  in  the  simulation.  If  one  adds  the 
approximate  energy  gained  from  the  well  to  the  maximum  kinetie  energy  assoeiated  with 
the  initial  Gaussian  wave  paeket,  an  approximate  value  of  the  maximum  momentum  ean 

be  determined  using  k^^^  » ^max  greater  than  the  eomputed 

value  the  momentum  wave  paeket  will  be  aeeurately  represented  throughout  the  entire 
simulation. 

The  parameter  At  is  important  in  the  error  assoeiated  with  the  Split-Operator 
method.  The  larger  the  time  step  the  less  time  required  to  propagate  the  wave  paekets  but 
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the  greater  the  error  assoeiated  with  the  split  operator.  As  noted  in  Appendix  A,  a  good 
rule  of  thumb  is  to  ensure  the  following  relationship  holds: 

E^At«h  (54) 

where  is  eomputed  by  adding  the  largest  kinetie  energy  assoeiated  with  the  initial 

Gaussian  wave  paeket  to  the  depth  of  the  well.  If  equation  (54)  is  maintained  the  error 
assoeiated  with  the  time  step  is  negligible  and  the  split-operator  will  remain  stable. 

The  propagation  times  ,  should  be  chosen  to  ensure  that  the  Moller  states  are 

max 

fully  formed,  and  that  the  time  dependent  correlation  functions  described  in  equation  (48) 
start  at  zero,  gradually  increase  as  the  Moller  states  overlap,  and  then  go  back  to  zero. 

The  energy  spacing  of  the  energy  grid  is  completely  characterized  by  the  momentum  and 
the  mass  of  the  particle  and  is  represented  as 

fi^kf 

EE  = - (55) 

2m 

The  spacing  of  the  energy  grid  is  not  uniform.  Although  a  Fourier  transform  is 
performed  to  convert  the  time  dependent  correlation  functions  into  energy 
representations,  the  requirement  that  the  spacing  be  uniform  is  only  required  for  a  FFT. 

It  is  important  to  understand  that  there  is  no  perfect  guide  to  choosing  parameters 
in  the  CPM.  The  relationships  described  above  are  meant  to  serve  as  a  guide  in 
formulating  the  algorithm.  An  example  of  the  parameters  used  to  perform  the  simulations 
in  this  thesis  is  provided  in  Appendix  C.  The  goal  is  to  define  parameters  which  decrease 
propagation  time,  minimize  grid  size,  and  accurately  represent  the  wave  functions 
throughout  the  simulation. 
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V  Properties  of  the  New  Algorithm 
5,1  Verifying  the  New  Approach 

Several  evaluations  were  performed  throughout  this  thesis  using  the  new 
application  of  the  CPM.  An  initial  simulation  was  performed  using  product  and  reactant 
states  having  only  positive  or  negative  momentum.  The  ri  -matrix  of  equation  (46)  is 
diagonal  in  this  case.  A  square  well  was  used  for  this  simulation  and  is  described  by 
equation  (31).  The  momentum  representations  of  the  initial  product  and  reactant  wave 
packets  were  centered  about  =  ±6 ,  similar  to  that  shown  in  Figure  3.  The  purpose  of 

this  simulation  was  to  prove  that  the  new  algorithm  produced  results  consistent  with 
previous  applications.  Figure  10  shows  the  results  of  the  correlation  functions  (real  part) 
defined  by  equation  (48).  Notice  that,  as  mentioned  previously,  the  correlation  functions 
start  at  zero  and  end  at  zero. 


Figure  10.  Real  part  of  the  correlation  functions  for  simulation  1 .  The  dashed  function  represents 
Cj(0=C3(0  ,  while  the  solid  line  represents  C2{t)  =  C^{t)  . 
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Figure  1 1 .  Transmission  (upper  )  and  reflection  (lower)  curves  for  the  square  well  of  depth 
V  =  —0.01  G.U  .  and  length  2ci  —  4  .  The  dashed  lines  represent  the  analytic  solutions 
while  the  solid  lines  indicate  the  numerical  results. 

As  shown  by  Figure  3,  the  wave  amplitudes  are  effectively  non  zero  for  the  momentum 
range  from  3  to  9  a.u.  Since  the  other  wave  packets  were  centered  about  =  -6  the 

corresponding  amplitudes  were  non  zero  for  the  momentum  range  from  -3  to  -9  a.u. 

Using  the  mass  of  the  particle,  1833.35  a.u.,  to  convert  the  momentum  range  to  an  energy 
range,  one  would  expect  good  results  from  approximately  0.0025  to  0.022  a.u. 

The  transmission  and  reflection  curves  are  shown  in  Figure  1 1  along  with  the 
analytic  solution  for  a  square  well.  Notice  that  the  results  begin  to  deviate  at 
approximately  the  energy  limits  given  above.  The  deviations  underscore  the  limitation 
that  only  S-matrix  elements  corresponding  to  the  range  of  momentum  for  which  the  wave 
packet  is  non  zero  are  obtainable.  After  verifying  the  algorithm  worked  correctly  we 
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began  investigating  scattering  models  for  which  the  product  and  reactant  wave  packets 
contained  a  significant  amplitude  for  both  positive  and  negative  momentum  ranges. 

5.2  Low-Energy  Scattering  (Square  Well) 

5,2.1  Difficulty  in  Constructing  the  Moller  States 

The  focus  of  the  thesis  is  to  investigate  scattering  when  both  positive  and  negative 
momentum  exist  in  the  product  and  reactant  states.  As  such  the  setups  depicted  in  Figures 
7  and  8  were  used  to  investigate  low  energy  scattering.  Figure  8  shows  a  small  overlap  of 
negative  and  positive  momentum  for  the  products  and  reactants.  Preliminary  calculations 
revealed  that  constructing  the  Moller  states  would  be  problematic.  The  trouble  is  that 
when  positive  and  negative  momentum  is  present  the  wave  packets  spread  in  both 
directions.  For  example,  consider  attempting  to  form  the  intermediate  Moller  state  for 
productl,  I  ,  shown  in  Figure  8  .  As  the  wave  packet  is  propagated  forwards  in  time 
under  most  of  it  will  move  to  the  right  leaving  the  interaction  region.  However,  a 

portion  will  also  move  to  the  left  since  it  also  contains  negative  momentum.  Movement  in 
both  directions  makes  it  very  difficult  to  have  the  wave  packets  completely  leave  the 
interaction  region  and  thus  accurately  form  the  Moller  states.  Consider  Figure  12. 

The  figure  shows  the  resulting  wave  packets  of  Figures  7  and  8  after  attempting  to 
construct  the  intermediate  Moller  states.  The  interaction  region  again  is  the  square  well 
described  in  equation  (3 1).  The  states  were  propagated  for  800000  atomic  units  and  yet 
portions  of  the  wave  packets  still  remain  within  the  square  well  region  (see  Figure  12b). 
Another  important  feature  of  Figure  (12)  is  the  size  of  the  coordinate  grid  needed  to 
represent  the  wave  packets. 
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Figure  12.  (a)  Approximate  Moller  states  after  t  —  +800000  (a.u.) 
propagation  time.  Notice  that  the  wave  functions  have  spread  dramatically  . 

(b)  Zoomed  in  portion  of  interaction  region  —2  <  X  <  2  . 

When  the  wave  paekets  are  eonstrueted  to  have  only  positive  or  negative 
momentum  a  grid  size  of  only  40  a.u.  is  needed  to  construct  the  intermediate  Moller 
states  (see  Figure  4).  However,  when  the  reactants  and  products  contain  an  overlap  of 
negative  and  positive  momentum  as  shown  in  Figure  8  the  grid  size  grows  to  10000  a.u. 
and  the  wave  packets  are  still  not  completely  outside  the  interaction  region.  According  to 
equation  (50)  the  larger  grid  then  requires  a  larger  N  value,  since  in  order  to  model  a 
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square  well,  Ax  must  be  very  small  (see  Figure  9).  The  larger  N  value  makes 
eomputation  diffieult  sinee  the  simulation  times  get  very  large.  Despite  these  drawbaeks, 
an  attempt  was  made  at  ealeulating  the  S-matrix  elements  using  the  approximate 
intermediate  Moller  states  shown  in  Figure  12. 

5.2,2  Approximate  Square  Well  Results 

The  simulation  was  performed  under  the  setup  deseribed  in  Figures  7  and  8  using 
the  approximate  intermediate  Moller  states  shown  in  Figure  12.  In  an  effort  to  have  the 
eomputation  time  reasonable  N  was  ehosen  to  be  2'®  for  this  simulation.  In  order  to 
represent  the  intermediate  Moller  states  of  Figure  12  a  grid  size  of  12000  atomie  units 
was  ehosen  leading  to  a  Ax  «  0.183  (see  equation  50).  The  large  eoordinate  grid  spaeing 
results  in  the  trapezoidal  potential  funetion  shown  in  Figure  13.  The  true  square  well  of 
equation  (31)  is  shown  in  the  figure  for  eomparative  purposes. 


X  {atomic  units} 


Figure  13.  Approximate  square  well  (dashed)  and  square  well  (solid) 
for 


Vix)  = 


J-O.Ol 

[o 


-2<x<2 

elsewhere 
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It  is  evident  that  the  results  will  eontain  error  since  the  potential  function  is  not  exactly  a 
square  well.  The  labeled  points  demonstrate  the  discreteness  of  the  numerical  potential. 
In  fact,  the  computer  is  representing  a  trapezoidal  potential  well  with  linear  regions 
joined  by  the  two  labeled  points.  However,  based  on  Figure  13  the  two  potential 
functions  are  reasonably  close  and  so  one  would  not  expect  to  see  an  extremely  bad  fit. 
The  results  of  the  simulation  are  shown  in  Figure  14.  Notice  that  the  simulated  results 
appear  shifted  from  the  analytic  solution  due  to  the  error  in  the  potential.  The  other 
distinct  feature  is  that  at  low  energies  the  simulated  results  begin  to  oscillate.  Based  on 
previous  calculations  it  seems  reasonable  to  assume  that  if  Ax  was  made  small  enough 
the  curves  would  line  up  and  thus  on  can  say  that  the  oscillations  appear  to  ring  about  the 
analytic  values.  One  should  pay  careful  attention  that  the  behavior  at  low  energies  here  is 
much  different  from  results  where  the  reactant  and  product  wave  packets  contained  only 
positive  or  negative  momentum.  The  simulated  results  are  oscillating  about  the  known 
values  rather  than  diverging  as  shown  in  Figure  1 1 . 

The  other  interesting  property  is  that  and  curves  oscillate  less  than 

theA^j.^^  and  curves.  The  reflection  coefficient  approaches  one  at  low  energies 

while  the  transmission  function  goes  to  zero.  The  empirical  observation  then  is  that  the 
amplitude  of  the  oscillations  near  low  energy  depends  on  the  magnitude  of  the  correct 
value  when  approximate  Moller  states  are  used.  In  order  to  perform  an  error  analysis 
showing  that  the  oscillations  diminish  and  the  results  converge  to  the  square  well 
solutions  one  would  propagate  the  reactant  and  products  for  longer  times  in  constructing 
the  Moller  states. 
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Probability  Amplitudes 


Figure  14.  (a)  Transmission  [  5^^  ,  5  ^  ^  ]  and  reflection!  5^^  ^ ,  5  ^  ]  curves  for  the  square  well  of 

depth  K  =  —0.01  a.U  .  and  length  2t3!  =  4  .  The  analytic  solution  s  are  (solid)  while  the  simulated  results 
are  (dashed),  (b)  Zoomed  in  portion  of  low  energy  region  showing  oscillatory  behavior. 


The  longer  the  reaetants  and  produets  are  propagated  the  further  they  move  outside  the 
interaetion  region  and  thus  the  closer  one  gets  in  obtaining  the  correct  Moller  states. 
Provided  Ax  is  made  small  enough,  a  plot  of  error  vs.  the  time  the  wave  packets  are 
propagated  away  from  the  interaction  region  should  reveal  convergence.  To  perform  the 
analysis  one  must  first  establish  a  measure  of  convergence  for  Ax  in  order  to  understand 
just  how  small  the  coordinate  grid  spacing  needs  to  be  to  approximate  the  analytic 
solution  of  the  square  well. 

5.2,3  Convergence  to  Square  Well  based  on  Ax 

In  order  to  establish  the  convergence  of  the  algorithm  four  tests  where  performed 
at  different  values  of  Ax  while  holding  all  other  parameters  constant.  The  relative  error 
was  computed  for  each  of  the  four  simulations  based  on  the  theoretical  transmission 
coefficient  for  a  square  well.  The  setup  depicted  in  Figures  2  and  3  was  used  to  perform 
these  tests.  A  plot  of  the  relative  error  from  a  sample  test  is  provided  in  Figure  15. 
Measurements  were  taken  at  a  “low”  energy  value  since  this  thesis  is  primarily  concerned 
with  low  energy  scattering. 


Figure  15.  Relative  error  vs.  Energy  for  sample  test  for  square  well  of 
depth  K(a:)=-0.01((3.w.)  and  width  la  =  A  .  The  energy  values  of  the 
measured  data  point  is  labeled. 
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Figure  16.  Relative  error  vs.  Ax  for  square  well  of  depth  and  width 

2a  =  4  .  An  approximate  fit  reveals  that  the  error  exponentially  increases  as  the 
coordinate  grid  spacing  is  increased. 

Notice  that  the  relative  error  approaches  small  values  as  Ax  is  deereased.  The  resultant  fit 
gives  one  an  idea  of  the  neeessary  Ax  required  to  aecurately  represent  the  square  well.  It 
should  be  noted  that  even  for  very  small  Ax  values  there  is  still  some  relative  error.  The 
explanation  for  this  is  that,  as  discussed  previously  there  is  error  assoeiated  with  At . 
Aeeording  to  Appendix  A,  the  error  term  associated  with  At  goes  like 


2 


(56) 


k  ^ 

where  =  Vo-\ — 2222^ .  For  this  ease  ~ll,  At  =  \,V^  =  0.01  and  m  =1833.35 , 

2m 

whieh  lead  to  an  error  term  of  approximately  0.001;  in  good  agreement  with  the  fit  of 
Figure  16  showing  that  the  error  goes  to  0.0019  as  Ax  ^  0 . 

According  to  Figure  16,  a  relative  error  of  approximately  0.5%  requires  a 
Ax  »  0.033  .  With  the  grid  size  used  in  the  diseussion  of  seetion  5.2.2 , 12000  a.u.,  the 
number  of  grid  points  needed  to  aehieve  a  Ax  »  0.033  would  be  2'‘^ .  When  the  wave 
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packets  contain  only  positive  or  negative  momentum  the  grid  size  is  only  40  a.u  (see 
Figure  4),  and  requires  N  =  2'^  to  achieve  the  same  grid  spacing.  Therefore,  compared  to 
considering  wave  packets  that  contain  only  positive  or  negative  momentum  to  those 
which  contain  both,  the  number  of  grid  points  needed  to  accurately  represent  the  square 
well  to  within  0.5%  relative  error  increases  by  a  factor  of  2^ .  As  mentioned  previously, 
the  parameter  N  is  important  in  considering  the  time  required  to  perform  a  simulation. 
The  number  of  computations  for  each  FFT  goes  as  N  log  N .  Therefore,  as  a  practical 
issue  it  is  beneficial  to  use  a  low  N  value  in  order  to  keep  the  simulation  times 
reasonable. 
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VI  Square  Well  Error  Analysis 

6,1  Methodology 

The  results  of  seetion  5.2.2  showed  that  an  error  analysis  of  the  observed 
“ringing”  was  needed  to  show  that  the  simulated  results  eonverged  to  the  eorreet  values. 
Using  the  parameters  N  =  2'*’ and  Ax  »  0.183 ,  we  deeided  to  eompare  several  tests,  eaeh 
eonsisting  of  larger  .  The  longer  the  produets  and  reaetants  are  propagated  into  the 

min 

asymptotie  ehannels  the  eloser  they  get  to  eompletely  leaving  the  well,  and  henee 
eonverging  to  the  eorreet  Moller  states.  The  problem,  however,  was  that  it  was  unelear  as 
to  what  the  “real”  solution  was  for  a  large  Ax  and  thus  there  was  no  baseline  solution  to 
eompare  the  results.  The  reader  should  be  reminded  of  the  diffieulty  in  formulating  a 
solution  when  Ax  is  large  (see  Figure  13).  At  large  Ax  the  diserete  representation  of  the 
potential  takes  on  a  trapezoidal  form.  An  estimated  relative  error  using  the  parameters 
N  =  2'®  and  Ax  «  0.183  was  ealeulated  to  be  approximately  37%  based  on  Figure  16. 
Therefore,  eomparing  the  results  to  the  analytie  square  well  of  equation  (31)  would  be 
eomputationally  prohibitive.  Having  a  baseline  solution  is  desired,  however  it  is  not 
neeessary  to  show  the  eonvergenee  of  the  eomputed  S-matrix  elements. 

Previous  plots  in  this  thesis  (such  as  Figure  11)  have  demonstrated  that  if  Ax  is 
made  small  enough  the  simulated  results  converge  to  the  correct  solution  up  to  a  certain 
point  at  very  low  energy.  Since  Figure  14  showed  that  the  results  “behaved”  like  the 
analytic  solution  across  the  entire  energy  range  it  was  no  longer  a  concern  to  make 

Ax  very  small.  Relaxing  the  requirement  of  small  Ax  allowed  me  to  explore  the  effects 
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of  increasing  without  having  to  have  unreasonably  long  simulation  times.  The 

min 

“ringing”  phenomenon  seen  in  Figure  14  is  not  dependent  on  Ax  but  eaused  primarily  by 
the  intermediate  Moller  states  not  clearing  the  interaetion  region  .  While  an  analytie 
solution  is  desired,  it  is  not  needed  to  show  that  the  amplitude  of  the  oseillations  at  low 
energy  converge  to  zero.  The  wave  paeket  residual  left  within  the  interaetion  region  is 
deereased  by  propagating  the  reaetant  and  produet  wave  paekets  for  longer  times. 
Therefore,  as  is  inereased  the  portions  of  the  wave  paeket  left  within  the  well 

min 

diminishes  and  the  eloser  the  states  beeome  to  the  eorreet  intermediate  Moller  states.  The 
assumption  made  in  eondueting  this  kind  of  error  analysis  is  that  the  results  are  eonsistent 
with  what  the  analytic  values  of  the  S-matrix  elements  would  be  provided  they  were 
known  for  the  ease  of  the  trapezoidal  potential  shown  in  Figure  13.  The  assumption, 
however,  is  not  unreasonable  eonsidering  that  previous  results  deseribed  in  this  thesis 
indieated  that  for  small  Av  the  eurves  behave  like  the  analytie  solutions.  The  results  of 
the  next  seetion  will  show  that  the  amplitude  of  these  oseillations  weakens  as  the  eorreet 
Moller  states  are  approaehed. 

6,2  Results 

The  time  for  whieh  the  states  were  propagated  was  inereased  from  t  =  4  x  10"^  to 
8x10^  atomie  units  in  inerements  ofS.OxlO"^  The  potential  used  to  perform  these 


^  Ringing  is  also  observed  when  the  eorrelation  funetion  is  not  allowed  to  go  baek  to  zero  and  at  energy 
values  not  defined  within  the  energy  spread  assoeiated  with  the  initial  wave  paekets. 
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simulations  is  depicted  in  Figure  13  as  the  dashed  trapezoidal  well.  As  discussed  earlier 
the  amplitude  of  the  ringing  should  diminish  the  longer  the  product  and  reactant  wave 
packets  are  propagated  away  from  the  interaction  region.  In  order  try  and  incorporate  any 
shift  in  phase  or  frequency  an  average  of  the  amplitude  of  five  oscillations  present  in  each 
case  was  measured.  An  example  of  the  oscillations  used  is  shown  in  Figure  17.  The 
general  trend  is  that  the  amplitudes  of  the  oscillations  are  smaller  at  higher  energies. 
Therefore,  if  one  can  show  that  the  average  of  these  oscillations  goes  to  zero  as  the  time 
of  propagation  is  increased  a  logical  assertion  would  be  that  the  oscillations 


Figure  17.  and  S_^.  curves  showing  the  oscillations  used  to  perform  the  error 

analysis. 


which  were  not  measured  decrease  as  well.  The  results  are  shown  in  Figure  18.  The 
average  amplitude  of  the  five  oscillations  was  measured  and  plotted  as  a  function  of  the 
time  for  which  the  products  and  reactants  were  propagated  during  the  formation  of  the 
approximate  Moller  states. 
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The  results  shown  in  Figure  18  underscore  the  assertion  that  the  oscillations 
decrease  the  further  in  time  the  product  and  reactant  states  are  propagated  away  from  the 
interaction  region.  A  linear  fit  to  the  data  gives  a  crude  indicator  of  the  time  needed  to 
diminish  the  ringing  to  some  reasonable  tolerance.  Assuming  that  the  simulated  results 
match  the  analytic  solution  for  the  trapezoidal  potential  well  used  to  perform  these 
calculations,  these  results  indicate  that  the  algorithm  converges  to  the  correct  S-matrix 
values  as  the  true  Moller  states  are  approached. 

It  would  then  be  reasonable  to  assume  that  if  the  trapezoidal  potential  well 
solution  was  known  the  computed  S-matrix  elements  would  be  close  to  the  analytic 
values.  However,  as  with  any  hypothesis  definitive  proof  is  necessary  to  validate  the 
assumptions.  As  such,  it  was  necessary  to  determine  the  analytic  transmission  and 
reflection  coefficients  for  a  trapezoidal  well. 


Figure  18.  Average  amplitude  of  the  five  oseillations  shown  in  Figure  17  vs.  the  time  for  whieh  the 
produet  and  reaetants  were  propagated  away  from  the  interaetion  region. 
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VII  Low-Energy  Scattering  (Trapezoidal  Well) 


7,1  Solving  Schrodinger’s  Equation 

Schrodinger’s  equation  for  a  trapezoidal  potential  well  eontains  the  potential, 


V  (x,  Ax)  = 


0 

-—(x  +  a) 
Ax^  ’ 

-K 

—  (x-a) 
Ax^  ^ 

0 


x<-a 

-a<x<-a  +  Ax 
-a  +  Ax<x<a-Ax 
a-  Ax  <  X  <  a 
x>a 


(57) 


and  is  shown  below. 


A 


-4 


I 


Figure  19.  Trapezoidal  potential  showing  effect  of  having  a 
large  Ax  .  The  five  regions  are  labeled  under  their  respective 
parts. 


The  standard  approaeh  is  to  mateh  boundary  eonditions  aeross  the  five  regions.  The 
solutions  in  regions  (I,  III,  and  V)  are  plane  waves  while  regions  (II  and  IV)  eontain  Airy 
funetions.  Mathematiea  was  used  to  symbolieally  solve  for  T  and  R  in  terms  of 
Ax,  a,  and  E  (Appendix  B).  The  transmission  eoeffieient  was  eheeked  by  eomparing  a 
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large  Ax  and  very  small  Ax  to  a  square  well  potential  of  length  2a  =  4 .  If  the  trapezoidal 
expression  is  eorrect  it  should  approach  a  square  well  solution  as  Ax  gets  small.  It  was 
also  verified  that  T  +  =  1 .  The  results  of  these  tests  are  shown  in  Figures  20  and  21. 

Comparing  Figures  20  and  21  it  is  apparent  that  the  trapezoidal  solution  approaches  the 
square  well  as  Ax  gets  small.  Figure  22  also  reveals  that  the  sum  of  the  transmission  and 
reflection  coefficients  is  equal  to  one.  These  observations  underscore  the  predicted 
expectations  of  the  correct  solution,  and  therefore  the  solution  was  taken  to  be  acceptable. 


(a) 


(b) 
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Figure  20.  (a)Plot  of  the  Transmission  coefficient  vs.  energy  (a.u.)  for  a  square  well  (solid 
line)  and  trapezoidal  well  (dashed)  with  Ax  =  0.1  (a.u.)  (b)  Plot  showing  the  comparison 
of  the  trapezoidal  potential  and  a  square  potential,  (c)  Relative  error  of  the  trapezoidal 
transmission  coefficient  and  square  well  vs.  energy  (a.u.) 
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(a)  (b) 


Figure  21.  (a)Plot  of  the  Transmission  coefficient  vs.  energy  (a.u.)  for  a  square  well 
(solid  line)  and  trapezoidal  well  (dashed)  with  Ax  =  S  (Machine  Epsilon),  (b)  Plot 
showing  the  comparison  of  the  trapezoidal  potential  and  a  square  potential,  (c)  Relative 
error  of  the  trapezoidal  transmission  coefficient  and  square  well  vs.  energy  (a.u.) 


T+R 
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Figure  22.  Plot  of  the  T+R  vs.  Energy  (a.u.)  showing  that  it 
indeed  equals  one  for  all  energies 
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VIII  Trapezoidal  Well  Error  Analysis 

Having  the  known  analytic  solution  for  a  trapezoidal  well  allows  a  large  grid  as 
well  as  a  large  Ax  value.  The  only  requirement  is  that  the  linear  regions  of  the  trapezoidal 
representation  are  sampled  well.  The  wave  packets  were  set  up  as  described  in  Figures  7 
and  8  except  with  the  trapezoidal  potential  shown  below 


X  {atomic  units} 

Figure  23.  Trapezoidal  potential  well  showing  discrete 
points  along  curve.  The  well  is  0.01  atomic  units  deep. 

Four  tests  were  run  at  increasing  times  to  determine  if  the  algorithm’s  results  converged 
to  the  analytic  solution  as  the  wave  packets  were  propagated  away  from  the  interaction 
region.  A  sample  calculation  showing  the  transmission  and  reflection  curves  is  shown  in 
Figure  24.  Notice  that  the  same  oscillatory  nature  seen  earlier  is  present.  However,  these 
results  are  different  than  those  of  Figure  14  in  that  the  simulated  results  are  nearly  an 
exact  match  outside  the  “ringing”  portions;  convergence  to  the  analytic  solution  would  be 
more  accurate  if  one  samples  more  points  along  the  linear  portions. 
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Figure  24.  (a)  Transmission  [  ]  and  reflection  [  ,  S  ]  curves  for  the  trapezoidal  well. 

The  analytic  solutions  are  (solid)  while  the  simulated  results  are  (dashed),  (b)  Zoomed  in  portion  of  low 
energy  region  showing  oscillatory  behavior. 


An  error  analysis  was  done  by  summing  the  differenee  between  the  simulated  results 
and  the  analytie  solution  over  the  region  for  whieh  ringing  was  observed: 
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(58) 


J=a  ''  ' 

and  >  (59) 

j=a  ''  '' 

where  £'^  =  1.15x10'^  and  =  2.5  x  1  O  '*  deseribe  the  first  and  last  elements  in  the  energy 

array  for  whieh  oseillations  were  measured,  R  and  T  denote  the  analytie  refleetion  and 
transmission  eoeffieients  for  a  trapezoidal  well,  and  and  dj  denote  the  summation  of 

the  error  assoeiated  with  refleetion  and  transmission.  The  portion  of  the  eurves  outside 
the  region  for  whieh  ringing  was  seen  was  not  ineluded  in  the  summation  sinee  the  error 
was  negligible  there. 

8,1  Results 

Measurements  were  taken  using  the  results  of  four  simulations  for  whieh  the 
produets  and  reaetants  were  propagated  away  from  the  interaetion  region.  A  table  of  the 
times  used  is  shown  in  table  2.  A  plot  of  5^^  and  dj  vs.  the  times  listed  in  table  2  is  shown 
in  Figure  25  along  with  a  trend  line  showing  that  the  summation  errors  go  as 

(60) 

Notiee  that  for  both  refleetion  and  transmission  the  error  eonverges  as  the  time  is 
inereased.  The  error  assoeiated  with  transmission  is  lower  than  the  refleetion  error  for  all 
times;  in  agreement  with  the  previous  error  analysis.  While  it  is  not  exaetly  elear  how  the 
error  made  in  the  Moller  states  eharaeterizes  the  error  made  in  the  transmission  and 
refleetion  eurves,  it  is  apparent  that  the  eloser  one  gets  to  making  the  eorreet  Moller  states 
the  eloser  both  eurves  beeome  to  the  analytie  solutions. 
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Test  1 

t  =  250000  a.u 

Test  2 

t  -  500000  a.u 

Tests 

t  =  750000  a.u 

Test4 

t  =  1000000  a.u 

Table  2.  Times  of  propagation  for  Error 
Analysis  Two 


t  {atomic  units} 


Figure  25.  (solid)  and  Sj,  (dashed)  vs.  time  along  with  power  fits  showing  convergence. 


To  test  the  assertion  that  the  ringing  seen  in  the  low  energy  portion  of  the 
transmission  and  refleetion  curves  is  caused  by  using  incorrect  Moller  states  we  measured 
the  area  of  the  residual  wave  packet  left  inside  the  interaction  region  at  the  time  the 
approximate  intermediate  Moller  states  were  formed.  At  each  time  listed  in  table  2  the 
area  under  the  curve  described  by  the  free  particle,  time-  evolving  Gaussian  wave  packet 
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described  by  equation  (26)was  measured  (Figure  26).  Since  the  analytic  form  was  known 
the  propagations  as  well  as  the  integrals  were  performed  analytically. 


Figure  26.  Schematic  showing  integrated  residue  (shaded)  inside 
interaction  potential  (Trapezoidal  Well). 


The  results  are  plotted  in  Figure  27  along  with  the  previous  plots  (5^  and  d-j. .  The 
integrated  data  was  scaled  such  that  the  point  at  t  =  1  x  1 0®  matched  (5^ (t  =  1x10'’)  and 
(5^(t  =  1 X  lO®) ;  the  same  scaling  factor  was  then  applied  to  the  other  three  data  points. 

Since  the  data  was  scaled  to  have  the  points  at  t  =  1  x  1  O'*  match,  the  integrated  data  is 
more  closely  related  to  the  relative  error  curves  at  longer  times.  If  the  points  at 
t  =  2.5  X 10^  were  scaled  to  match  the  plots  would  indicate  a  stronger  correlation  at 
shorter  times.  Although  the  scaling  is  arbitrary  the  empirical  results  shown  in  Figure  27 
underscore  the  assertion  that  the  error  seen  in  the  ringing  is  related  to  the  residual  tail  left 
inside  the  interaction  region  during  the  construction  of  the  Moller  states. 
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Figure  27.  5^  &  5j.  along  with  the  integrated  data  showing  the 

correlation  between  residue  left  inside  interaction  region  and  error 
associated  with  ringing. 

The  most  probable  explanation  for  the  offsets  is  that  the  error  associated  with  the 
incorrect  Moller  states  is  propagated  through  the  rest  of  the  CPM  method  (ie.  correlation 
function,  Fourier  transforms,  and  the  S-matrix  elements  ).  Both  curves  show  the 
drawback  of  this  approach  in  that  it  requires  long  propagation  times  to  achieve  little  error 
in  the  low  energy  regions. 
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IX  Conclusion 


This  thesis  contributes  to  the  understanding  of  the  Channel  Packet  method  and 
low-energy,  time  dependent  scattering  theory.  Two  properties  of  this  thesis  make  it  ideal 
for  understanding  the  complications  associated  with  the  low  energy  problem.  The  first  is 
that  it  demonstrates  the  need  for  a  four-by-four  matrix  of  expansion  coefficients  to  handle 
wave  packets  which  contain  positive  and  negative  momentum.  In  doing  so  however,  one 
must  accept  the  increased  complexity  of  solving  all  four  S-matrix  elements  at  once.  The 
second  property  is  that  this  thesis  explains  the  difficulty  in  constructing  the  Moller  states 
from  reactant  and  product  wave  packets  which  contain  both  positive  and  negative 
momentum. 

The  thesis  extended  the  application  of  the  CPM  to  tackle  the  specific  problem  of 
low  energy  scattering.  A  new  algorithm  was  presented  which  relaxes  the  requirement  that 
the  product  and  reactant  wave  packets  only  contain  positive  or  negative  momentum.  In 
doing  so,  the  algorithm  allows  wave  packets  to  be  fully  represented  across  the  entire 
energy  range.  The  requirements  to  solve  the  low  energy  problem  are  now  completely 
understood.  The  first  requirement  is  that  the  product  and  reactant  wave  packets  contain 
both  positive  and  negative  momentum.  The  second  requirement  is  that  the  correct  Moller 
states  must  be  used.  While  the  technique  presented  in  this  paper  meets  the  first 
requirement,  the  latter  was  never  obtained.  However,  this  thesis  demonstrated  that  as  the 
proper  Moller  states  are  approached  the  oscillations  diminish  and  the  correct  S-matrix 
elements  are  approached.  It  was  shown  that  to  accurately  construct  the  Moller  states  will 
require  a  very  large  grid,  and  that  the  main  downside  to  the  approach  taken  in  this  thesis 
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is  that  one  must  accept  very  long  simulation  times.  However,  as  technology  increases  and 
processors  become  faster  this  algorithm  may  become  the  algorithm  of  choice  for  time 
dependent  scattering  theory.  This  thesis  explained  the  difficulties  associated  with  low 
energy,  time  dependent  scattering  theory  as  well  as  a  new  approach  which  significantly 
improves  the  calculation  of  the  low-energy  S-matrix  elements.  The  thesis  demonstrated 
that  the  new  approach  will  work  but  at  a  high  computational  price.  Finally,  if  one  is  only 
interested  in  high  energy  scattering  using  the  previous  approaches  described  in  [1-3]  is 
simpler.  However,  if  one  needs  to  characterize  low  energy  scattering  the  algorithm 
described  in  this  thesis  is  a  valid  approach. 

Ideas  for  further  research  include  investigating  how  one  could  avoid  using  very 
large  grids  when  the  products  and  reactants  contain  positive  and  negative  momentum. 

One  possible  solution  may  be  that  during  the  second  half  of  the  Moller  operator,  the 
propagations  under  the  full  Hamiltonian,  the  Gaussian  wave  packets  could  be  “created”  at 
the  boundaries  of  a  smaller  grid  at  each  time  step  as  they  moved  back  into  the  interaction 
region.  The  split  operator  would  need  to  be  modified  to  accept  certain  boundary 
conditions  at  the  ends  of  the  grid.  If  one  could  construct  the  algorithm  using  small  grids 
the  computation  time  needed  to  perform  the  simulations  would  be  greatly  reduced. 
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Appendix  A.  The  Trotter  Theorem  (Split-Operator  Approximation) 

The  time  evolution  operator 
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0 .  In  order  to  split  the  operator  one  must  approximate  a  new  operator 


using  the  Trotter  theorem,  whieh  states  that  [13] 
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where  T,  and  are  general  Liouville  operators.  If  one  allows  M  to  be  large,  but  finite  : 
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whieh  ean  be  rearranged  by  taking  the  root  to  give 
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Interpreting  t/M  as  a  single  time  step.  At ,  and  replaeing  T,  and  Tj  with 


V  f 

- and - a 

n  n 


single-time-step,  approximate  propagator  may  be  written  as 
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This  operator  is  unitary  and  thus  preserves  time  reversibility: 

P\At)  =  P{-At)  =  P  \At)  , 

Pi-At)PiAt)  =  I. 

Using  a  Taylor  expansion  P{At)  is  aeeurate  to  order  Ap : 
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Henee,  P{At)  =  U{At)  up  to  seeond  order  in  At .  One  can  then  interpret  that  the 
propagator  is  stable  under  the  requirement 

,  (A.10) 

where  =  T  _  +  U  _  . 
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Appendix  B,  Trapezoidal  Well  Transmission  and  Reflection  Coefficients 
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Appendix  C  (Simulation  Parameters) 


•^min  -+1200 
max 

m  =  1833.35 

=  1 

Ax  «  0.092 

=  0.5 

A^  «  0.00026 

F„=-0.01 
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-  0.03 

1.12x10" 

Table.  3  Sample  parameters  of  simulation  for 

K=±A 
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